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ABSTRACT 

We propose a method to generate ‘genetically-modified’ (GM) initial conditions for high- 
resolution simulations of galaxy formation in a cosmological context. Building on the 
Hoffman-Ribak algorithm, we start from a reference simulation with fully random initial con¬ 
ditions, then make controlled changes to specific properties of a single halo (such as its mass 
and merger history). The algorithm demonstrably makes minimal changes to other properties 
of the halo and its environment, allowing us to isolate the impact of a given modification. As 
a significant improvement over previous work, we are able to calculate the abundance of the 
resulting objects relative to the reference simulation. 

Our approach can be applied to a wide range of cosmic structures and epochs; here we 
study two problems as a proof-of-concept. First, we investigate the change in density profile 
and concentration as the collapse time of three individual halos are varied at fixed final mass, 
showing good agreement with previous statistical studies using large simulation suites. Sec¬ 
ond, we modify the z = 0 mass of halos to show that our theoretical abundance calculations 
correctly recover the halo mass function. The results demonstrate that the technique is robust, 
opening the way to controlled experiments in galaxy formation using hydrodynamic zoom 
simulations. 
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1 INTRODUCTION 

Understanding galaxy formation requires us to take account of the 
variety of halo assembly histories that build today’s population. 
Many pressing questions — such as the origin of varying mor¬ 
phologies (e.g. van Dokkum et al. 2013; Papovich et al. 2015) and 
bulge sizes (Kormendy 2015) — will be answered by understand¬ 
ing the interplay between complex, non-linear physics and the var¬ 
ious histories for mass accretion. The fundamental difficulty is that 
these histories are in turn determined by the random initial condi¬ 
tions seeded in the early universe. 

This paper is the first in a series to directly tackle that prob¬ 
lem using a novel approach. The most typical solution is to sim¬ 
ulate large numbers of galaxies in a representative volume (e.g. 
Genel et al. 2014; Schaye et al. 2015; Codis et al. 2015). However, 
this is computationally expensive and limits the resolution that can 
be achieved for any single object. Conversely, zoom-in simulations 
achieve the maximum level of physical detail for a given compu¬ 
tational time. They have been used to establish that qualitatively 
different processes come into play at sub-kpc resolutions, where 
processes within the interstellar medium begin to be resolved (e.g. 
Governato et al. 2007; Guedes et al. 2011; Brook et al. 2011; Hop¬ 
kins et al. 2013; Pontzen & Governato 2014). But such approaches 
only sample over a small and potentially biased range of merger 
histories. A third tactic is to use isolated, idealised set-ups to test 


particular hypotheses (e.g. Naab et al. 1999; Robertson et al. 2006; 
Hopkins et al. 2009); but these by definition lack a full cosmolog¬ 
ical environment. Thus, it is difficult to quantitatively connect the 
results of isolated and zoom simulations to the observed galaxy 
population. 

Our aim is to combine the best aspects of these three types of 
numerical study. We proceed by systematically changing aspects 
of individual galaxies’ histories (such as mass and merger history) 
within a cosmological simulation, while keeping track of the sta¬ 
tistical likelihood of the changes to understand the relative abun¬ 
dance of objects of different types. This can be achieved by using 
the Hoffman-Ribak algorithm (Hoffman & Ribak 1991, hereafter 
HR91; see also Bardeen et al. 1986; Bertschinger 1987; van de 
Weygaert & Bertschinger 1996 for further theoretical background). 
A more common use for HR91 is to obtain simulations resembling 
the local universe by turning a given observational dataset (e. g. 
the local distribution of galaxies) into a prescription for the initial 
conditions of a numerical simulation (Bistolas & Hoffman 1998; 
Mathis et al. 2002; Kravtsov et al. 2002; Klypin et al. 2003; HeB 
et al. 2013; Jasche & Wandelt 2013; Sorce et al. 2014). There is a 
significant literature that uses this technique to study the formation 
history of the Local Group (Zavala et al. 2009; Klimentowski et al. 
2010; Libeskind et al. 2010; Iliev et al. 2011; Forero-Romero et al. 
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2011; Kitaura 2013; Doumler et al. 2013; Dayal et al. 2013; Nuza 
et al. 2014; Brook et al. 2014). 

Instead, we propose to use constrained initial conditions as an 
experimentation toolkit for the formation of a particular halo em¬ 
bedded in a cosmological volume. This approach has precedent: 
for example, Frenk et al. (1999) used the HR91 method to create 
galaxy cluster initial conditions for a comparative study of numer¬ 
ical simulation codes. More recently, Romano-Diaz et al. (2006, 
2007) and Hoffman et al. (2007) simulated a single dark matter ob¬ 
ject of ~ 10^^M© with different substructures to understand the 
impact of quiescent and violent accretion phases on the inner prop¬ 
erties of the halo, and the origin of the universal halo density pro¬ 
file. By including baryonic physics, Romano-Diaz et al. (2011a,b, 
2014) studied galactic properties in extremely overdense regions 
which may host the early precursors of QSOs. In a similar vein, 
Dubois et al. (2012) use the numerical implementation from Prunet 
et al. (2008) to investigate the accretion of material in the cores of 
very massive halos to shed light on the formation of black holes at 
high redshifts. 

In all of the above cases, a simulated object was created by 
constraining the properties of a region defined by an analytical pro¬ 
file (typically a Gaussian peak). Constrained properties included 
the height of the density peak at the origin and its first- and second- 
order derivatives (see e. g. Prunet et al. 2008 and the Appendix 
of Romano-Diaz et al. 2011a). This creates objects that are well- 
defined in a theoretical sense (e. g. one can predict their collapse 
time reasonably well), but that represent configurations which may 
or may not be common in fully random initial conditions (ICs). 

What sets our work apart from these previous efforts is that 
we always start with a ‘reference’ halo from a simulation based on 
fully random ICs. We are able to impose constraints on volumes 
of completely arbitrary shape, using the particles that make up a 
single dark matter halo embedded in a cosmological volume. Once 
the constraints are applied, we re-run the simulation and compare 
the results to the original reference run. 

This has two immediate benefits. First, we can fine-tune se¬ 
lected properties of the halo while demonstrably ensuring that the 
constrained object is as similar as possible to the reference run — a 
controlled ‘genetic modification’ (GM) of the halo. Second, we can 
calculate the change in the likelihood of the field after the modifi¬ 
cation; in other words, we can assess the relative abundance of the 
genetically-modified systems compared to the original. This will 
allow us to test whether connections between merger history and 
morphology quantitatively account for observed population statis¬ 
tics. 

The current work provides a first illustration of both these as¬ 
pects of the technique. Specifically, we study the properties of sev¬ 
eral halos as their total mass and merger history are systematically 
changed. We investigate the concentration at z = 0 for different 
mass accretion histories and find overall excellent agreement of 
our constrained halos with relations derived from statistical aver¬ 
ages over large simulations. There are many studies that connect 
the concentration parameter to other halo properties like the mass, 
collapse time or mass accretion history, halo shape and angular mo¬ 
mentum (e. g. Bullock et al. 2001; van den Bosch 2002; Wechsler 
et al. 2002; Zhao et al. 2003; Reed et al. 2005; Bett et al. 2007; 
Maccio et al. 2007; Neto et al. 2007; Duffy et al. 2008; Maccio et al. 
2008; Zhao et al. 2009; Ragone-Figueroa et al. 2010; Prada et al. 
2012; Ludlow et al. 2013, 2014; Klypin et al. 2014; Correa et al. 
2014,2015a,b). Often, these studies operate by considering a statis¬ 
tical sample from a large volume simulation to find correlations and 
provide fitting functions. Even though the statistical power in recent 


simulations is excellent, there is still considerable scatter around the 
median relations. Since the density profile of dark matter halos is 
an important ingredient in theoretical models, it is important to un¬ 
derstand these correlations and the scatter. Given the large number 
of parameters that could infiuence the evolution of a halo, princi¬ 
pal component analysis has been used to investigate correlations 
between them (Skibba & Maccio 2011; Jeeson-Daniel et al. 2011; 
Wong & Taylor 2012). Our approach of designing ‘experiments’ in 
galaxy formation provides a complementary approach to computa¬ 
tionally expensive statistical studies. 

This paper is organised as follows: in Sec. 2 we give a brief 
outline of the HR91 technique and our specific implementation. 
Section 3 contains details of the numerical simulations that are used 
to obtain the results in the rest of the paper. In Sec. 4 we provide 
a brief illustration of some of the constraints we have applied to 
the reference initial density field, focusing on infiuencing a single 
halo traced by its particles. Next, we study the results of designing 
different merger histories for a set of halos in Sec. 5, focusing on 
their collapse-concentration relation. In Sec. 6 we discuss a method 
for assessing the relative abundance of the modified halos by defin¬ 
ing a measure, and show that our results are consistent with the 
cosmological halo mass function. We summarise in Sec. 7. Finally, 
Appendices A and B contain the mathematical details of our re¬ 
formulation of the HR91 technique including a translation between 
our notation and theirs. 


2 OUTLINE OF THE METHOD 

We now present a brief outline of the mathematical technique by 
which initial conditions can be generated that satisfy certain con¬ 
straints, while remaining consistent with a ACDM power spectrum. 
This technique is described in a slightly different formulation by 
HR91. The full derivation can be found in Appendix A. 

By assumption, the density field in the early universe is lin¬ 
early perturbed around a background density po, so that 

p{x) = po{l +5{x)), (1) 

where 5{x) is a Gaussian random field with statistical properties 
specified by the ACDM transfer function and inflationary tilt. 

Generating ICs involves sampling the Gaussian random field 
at a list of discrete points Xj,, where the integer value v de¬ 
cides which point we are discussing. In particular, when running 
a uniform-resolution cosmological simulation with N particles on 
a box side, v runs from 1 to . The sampled field then consists of 
an -length vector 5, where an element is given by di, = 5{xjp). 
The values of 5 are drawn from a multivariate Gaussian proba¬ 
bility distribution with mean {6) — /liq and covariance matrix 
Co = {{5 — — fjLo)). Cosmological initial conditions have 

zero mean, /lIq = 0, but the HR91 technique is not limited to this 
case. 

A constrained field is defined by requiring = (7 for some 
constraint vector a, which also contains elements. In general, 
d is real, though the formalism also extends to the case where it 
is complex-valued. A simple example would be to fix the density 
contrast to 0 at position xi. This requires = 1 for z/ = 1 (before 
normalisation, see below) and 0 otherwise, and (7 = 0. Throughout, 
we will use Greek indices in this way to denote values of either 5 
or a at a specific grid position Xu. 

To actually create a field satisfying any given constraint, one 
could sample repeatedly from the underlying population until ob¬ 
taining a realization that satisfies (or is close to satisfying) the 
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Name 

Halo 

Constraint 

Section 

Reference 

N/A 

none 

3 

H24-MH 

Halo 24 

^ 10 /^ref — { 0 . 5 , 1 . 5 } 

4,5 

H24-MH* 

Halo 24 

d{z = l)Kef = {0.9, 1.1} 

4,5 

H37-MH 

Halo 37 

^lo/^ref — {0.5, 1.5} 

4,5 

H37-MH* 

Halo 37 

d(z = l)Kef = {0.9, 1.1} 

4,5 

H40-MH 

Halo 40 

rfio/dref = {0.5, 1.5, 2} 

4,5 

H40-MH* 

Halo 40 

d{z = l)/dref = {0.9,1.1,1.2} 

4,5 

H24-mass 

Halo 24 

d/d^ei = {0.5,0.8,1.2,1.5} 

4,6 

H37-mass 

Halo 37 

d/dref = {0.5,0.8, 1.2, 1.5} 

4,6 

H40-mass 

Halo 40 

d/d^^{ = {0.5, 0.8,1.2,1.5} 

4,6 


Table 1. Overview of the simulations used in this paper. For more details on 
the individual runs, see the text in the sections mentioned in the last column. 
MH and MH* stand for the two different ways of constraining the merger 
history of the halo. The simulations marked in bold are not actually used 
because they are not in equilibrium at 2 ; = 0 (see Sec. 5 for a discussion). 

requirement. However, such an accept-reject algorithm would be 
computationally expensive to implement in practice; instead, the 
HR91 technique makes a mathematical rearrangement that requires 
only one set of random numbers to be generated. As detailed in Ap¬ 
pendix A, this rearrangement also shows that a constrained Gaus¬ 
sian random field remains Gaussian. This allows us to apply a large 
number of constraints independently, with the final result obeying 
ol\ 5 — di for each i where the Roman index i denotes the n dif¬ 
ferent constraints. The properties of the constrained field are then 
determined by a new mean and covariance 

n 

Mn = Mo+ X^Coai (di - aJ/Lig) (2) 

i=l 

n 

Cn = Co — ^ Coa^aJCo, (3) 

i = l 

provided that the {az} have been orthonormalized^ in the sense 
that cxlCocxj = 6ij. 

A realization of the constrained field could therefore be ob¬ 
tained by calculating and Cn and drawing random numbers 
accordingly. However, in practice, dealing directly with Cn from 
Eq. (3) becomes prohibitively expensive for large N. The problem 
is that, whereas Co is the ACDM power spectrum and therefore di¬ 
agonal in Fourier space, Cn is generally not sparse in either pixel 
or Fourier space. Instead, one can make the ansatz that a realization 
obeying n constraints, Sn, can be obtained starting from a realiza¬ 
tion of the unconstrained field, ^o, via 

dn — {do — + ^6^, (4) 

where Pn is a matrix that depends on Co and {a^}. 

By requiring that ^n obeys the correct statistics and, addition¬ 
ally, requiring that the changes made to the field are minimal, one 
can uniquely derive the HR91 solution for Pn - The details are given 
in Appendix A, with the result that 

n 

Sn — do ^ ^ doOli {di d-io) 5 (5) 

i = l 

where we have defined dio = aj^o to represent the value of the 
constrained quantities in the unconstrained realization, and again 

^ Note that this orthonormalization can always be arranged for any set of 
non-conflicting original constraints. 


require that the {az} are orthonormalized. This reduces the actual 
calculation to a series of vector multiplications and summations in 
Fourier-space (since Co is diagonal there). The memory require¬ 
ments are managable since we need only store vectors of length 
iV^, instead of the x matrix Cn- 

Consequently, a continuum of constrained realizations can be 
generated from a single realization of the original ensemble. We 
select a single dark matter halo from a reference run at z = 0, 
and then return to the initial conditions and place constraints on 
the particles that make up this object. In this way our constraint 
regions are defined directly via the halo particles, without requiring 
any assumptions about the properties of density peaks in the initial 
conditions, or any type of direct smoothing of the density field (only 
indirectly through the halo finding at 2 ; = 0). 


3 SIMULATION SETUP 

All our simulations were run with P-Gadget-3 (Springel 2005; 
Springel et al. 2008). The initial conditions have been set up at red- 
shift z = 99 and evolved to z = 0, saving 100 snapshots from 
z = 9 equally spaced in scale factor. The cosmological model is 
Qm = 0.279, Qb = 0.045, Qa = 0.721, erg = 0.817, h = 0.701, 
Us = 0.96, corresponding to a WMAP5 cosmology (Dunkley et al. 
2009). While these cosmological parameters have been revised in 
more recent datasets, they allow an easier comparison of our results 
with the literature. All simulations have a (comoving) box size of 
L = 50 h~^ Mpe ~ 71.3 Mpe, and iVpart = 256^ dark-matter 
particles, resulting in a particle mass of 8.24 • lO^M©. The Plum¬ 
mer equivalent force softening length (which limits the smallest 
accessible scales) is e = 25.6 kpc in comoving units. 

We use the SUBFIND code (Springel et al. 2001), which finds 
halos with the friends-of-friends (FoF) method. SUB FIND also iden¬ 
tifies subhalos inside the top level FoF groups, but we always use 
the whole group for particle tracking. Each FoF group is assigned a 
unique number, sorted in descending order by mass. SUB FIND also 
provides a list of halo particle IDs, which allows us to track the 
halos between snapshots and across different simulations (by se¬ 
lecting those objects which have the most particles in common at 
^ = 0). We choose the standard FoF linking length of 0.2 times the 
mean interparticle distance. 

Our analysis makes use of the Python module PYNBODY 
(Pontzen et al. 2013). We select halos with mass M 200 ~ lO^^M© 
at 2 : = 0, which are well-resolved but not the most massive (and 
therefore rare) objects in the box. Throughout this paper, M 200 
refers to the halo mass contained in r 2 oo, the radius within which 
the mass density is 200 times the critical density of the universe at 
that time. We construct halo merger trees by tracing halo particles 
from z = 0 backwards in time through each simulation snapshot. 
This allows us to determine the mass accretion history and other 
internal properties of each object as a function of time. We take ad¬ 
vantage of the fact that the density field is first set up by assigning 
one particle to each grid node (the displacements are applied later). 
This means that any particle at z = 0 can be traced back to a grid 
position Xjy in the initial conditions, and no additional interpolation 
is necessary to calculate S{xu). 

Table 1 gives an overview of the different runs that are used 
in this paper; the last column provides the section where the con¬ 
straints are described and the results are discussed. In order to show 
that the technique is robust, we have selected three halos of similar 
mass in the reference run, and constrained their properties in differ- 
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Figure 1. Left panel: The density of the reference ICs (black circles) and modified H40-MH-2 ICs (red crosses) for the early collapse constraint, where 
the density of the 10% innermost particles is increased by a factor of 2. The slice is 5 kpc wide in the y- and ^-coordinates, to give an impression of the 
3D structure. Each symbol corresponds to a single particle/initial grid point. The constrained density field maintains the complicated (sub-)structure that was 
present in the reference run. Right panels: The same two ICs as a 2D projection in the x — y-plane. Only those particles that form each halo at 2 ; = 0 are 
shown here; it is these particles that are used for generating the constraint in our algorithm. The higher central density is clearly visible in the constrained case. 
The results of these simulations will be discussed in detail in Sec. 5. 


ent ways. Therefore each simulation name contains a halo number 
and constraint type. 


4 ILLUSTRATION OF CONSTRAINTS 

We now present a simple illustration of the technique with which 
we generate a density constraint. We will discuss the results ob¬ 
tained from running simulations with these constraints in the fol¬ 
lowing sections. 

Our approach constrains the actual Lagrangian region that 
collapses into a halo at z = 0; by contrast, in previous work 
the constraints were typically chosen to follow some analytical 
form, in order to connect the constraints to theoretical models such 
as Press-Schechter theory (e. g. van de Weygaert & Bertschinger 
1996, Romano-Dfaz et al. 2006). Since we know which particles 
are going to collapse in the reference run, we do not need to as¬ 
sume a specific form for the peak or a smoothing scale. The result¬ 
ing constrained halo will be very similar to the reference object, 
unless the constraint radically changes the collapsing region, e. g. 
by introducing a large overdensity in a region which only forms an 
intermediate mass halo in the reference run. 

Designing the constraints for a given modification to the final 
halo requires a physical understanding of the evolution. Ultimately 
a proposed constraint must be tested by trialling the changes and 
testing that they have the desired effect and that they are statistically 
consistent with the modified halo existing in the unconstrained uni¬ 
verse. We will demonstrate both of these properties over the re¬ 
mainder of the paper. 

Changing the mass can be achieved by changing the density 
contrast of the halo particles in the initial conditions. By creating 
a larger or smaller overdensity, we influence the final mass by in¬ 


creasing or decreasing the overall size of the region which has the 
average threshold density to collapse by a specified redshift (Press 
& Schechter 1974). We term this a density constraint: in the initial 
conditions, we calculate the average mass overdensity of all A^part 
particles in the reference halo 


1 

N part 


-^part 

yf 5{x^) = d, 

P=1 


( 6 ) 


where we again use the fact that each particle corresponds to a grid 
position Xu. Before orthonormalization, the value of the constraint 
vector a is then 1/A^part for each particle which belongs to the 
halo, and 0 otherwise. The density can now be increased or de¬ 
creased by enforcing the value of d. Results of simulations with 
different choices for d that produce halos with higher or lower mass 
at z = 0 will be used in Sec. 6. 

More specifically, according to the Press-Schechter argument, 
the collapse time of a halo is related to its peak height 0 = 6 /cr(R), 
where (j{R) is the variance of the density field smoothed on a scale 
R. Therefore by fine-tuning the overdensity on different scales 
within the initial conditions we can modify the accretion history. 
In particular, by increasing (or decreasing) the density contrast in 
an inner region of the halo, while requiring that the overall density 
contrast of the halo particles stays the same, we are able to gener¬ 
ate a halo with very similar mass ai z — 0, but a faster (or slower) 
accretion history. 

Figure 1 illustrates an example of such a constraint acting on 
the initial conditions. In the left panel, we show a slice through the 
dark matter density field in the initial conditions, centred on the 
halo’s centre of mass. Each black circle corresponds to a density 
value in the reference run, and the red crosses show the same po¬ 
sition in the constrained run. The slice is 5 kpc wide in the y- and 
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Figure 2. Left panel: Mass accretion history for early (red solid) and late collapse (blue dashed) runs, expressed by the FoF mass (all particles assigned by 
the halo finder). The black solid line with points shows the same halo in the reference run; each point is one snapshot, illustrating the time resolution of our 
simulations. Right panel: same but for the virial mass M 200 ? which does not converge to a common value at late times because M 200 probes the inner regions 
of the halo (see Figure 3), which are affected by the collapse time. 


^-coordinates, to give an impression of the 3D structure. Here, we 
show the constraint that will be used in the ‘early’ run (discussed 
in detail in Sec. 5), where we have increased the density of the in¬ 
nermost 10% of halo particles to be a factor of two higher, while 
keeping the overall density the same. An alternative approach is to 
actually identify substructures at an intermediate redshift and apply 
the ‘inner’ constraints to those specific particles. We have tried this 
second approach for the present work, selecting the particles which 
have already collapsed around z ~ 1 by constructing a merger tree 
from the SUB FIND output in the reference run. Table 1 contains an 
overview of all simulations used in this paper; we denote the first 
method of constraining the inner region with diQ and the second by 
d{z = 1) there. The two modification methods give near-identical 
results (Sec. 5), making the outcomes reassuringly insensitive to 
the intuition guiding the modifications. We have found that it is also 
possible to add further constraints to modify the build-up in differ¬ 
ent subhalos and so fine-tune the accretion history to any required 
degree. 

As explained in Sec. 2, the modified field is constructed to fol¬ 
low the peaks and troughs of the underlying density field, thereby 
maintaining the same substructure as much as possible. In the right 
panel, we show the 2D projection (x — t/-plane) of the density of 
halo particles in the initial conditions, again for both the reference 
run and constrained run. In both cases, the density is calculated for 
all particles that are part of the halo at z = 0, which have been 
traced back to z = 99. The effect of increasing the density in the 
innermost region can be clearly seen in the constrained run (red 
border). In addition the second constraint, which keeps the overall 
mass the same, leads to a compensation effect in the initial condi¬ 
tions, removing some particles in the outer regions which fall into 
the reference halo but not the constrained one. We will discuss the 
results of simulations with these modified initial conditions in the 
next section. 


5 MERGER HISTORY AND 

CONCENTRATION-COLLAPSE RELATION 

So far we have looked at how applying various constraints mod¬ 
ifies the initial linear overdensity field. We will now consider the 
changes that result when the new initial conditions are used in a 
numerical simulation, starting with our modified merger history. 

The mass accretion histories for simulations with the H40- 
MH-2 ‘early’ (circles) and H40-MH-0.5 ‘late’ collapse constraint 
(crosses) are shown in Fig. 2. Here, we chose the innermost 10% 
of particles and changed their overdensity by a factor of 2 (0.5) for 
the early (late) collapse cases. In the left panel we show the time 
evolution of the total mass (including all substructures) for the two 
constrained runs and the reference halo. It is clear that the accretion 
rates differ quite significantly at early times, but are compensated 
at late times, leaving the overall mass of the objects the same. In the 
right panel we show the time evolution of M 200 . This quantity only 
measures the mass up to r 2 oo instead of the total mass of linked 
FoF particles (which can extend out to several r 2 oo). As in the for¬ 
mer case, the accretion rate follows the expected behaviour in the 
early and late collapse cases. However, at late times, M 200 differs 
between the constrained runs and the reference runs. This is due to 
a change in the halo density profile related to the collapse time, as 
we will now explore. 

The halo radial density profiles for these three simulations at 
z = 0 are illustrated in Fig. 3, with inset panels showing projected 
density maps. As the collapse is delayed, the slope in the inner re¬ 
gions becomes less steep. The location of the virial radius r 2 oo is 
indicated by an arrow in each case; as expected given our discus¬ 
sion above, this is displaced inwards by the relative shallowness of 
the late collapse H40-MH-0.5 case. 

The difference between the density profiles can be encapsu- 
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Figure 3. Density profile of the reference halo (black dot-dashed) and the 
‘early’ (blue dashed) and ‘late’ (red solid) constrained runs at 2 ; = 0. The 
leftmost arrow indicates the softening length of the simulation, and the other 
arrows indicate the virial radius of each halo. Inset panels: density projec¬ 
tion {x — y-plane) of the resulting halos at 2 ; = 0. All panels show a region 
2.5 Mpc across, include only the FoF group particles, and use the same 
colour scale for the column density. 


lated in the concentration parameter 

7^200 



(7) 


where r 2 oo is the virial radius of the halo (defined in Sec. 3), and rg 
is the scale radius in the NFW density profile (Navarro et al. 1997) 


p(r) = 


4ps 

(r/rs)(l + r/rs)2’ 


( 8 ) 


We fit an NFW profile to each of our halos at 2 ; = 0, after determin¬ 
ing its centre using a shrinking-sphere method and estimating the 
density in 100 radial bins of equal size. In order not to contaminate 
the fits with numerical artefacts introduced by the finite particle 
resolution, we exclude from the fit the innermost regions (2 times 
the softening length, e) which are affected by the force softening 
(e. g. Power et al. 2003), and regions with r > 0.6 r 2 oo, which 
may not be relaxed. We tested that the choice of the minimum and 
maximum radius has negligible impact on the estimate of rg. The 
measured values of the concentration are 4.2, 6.1 and 11.1 for late, 
reference and early simulations respectively. 

The GM method allows us to study the relationship between 
the collapse time of a halo (defined below) and its concentration as 
measured at z = 0. Using a large statistical sample, Wechsler et al. 
(2002) found that c oc scale factor at collapse time. We 

follow their procedure to obtain the collapse scale factor by fitting 
the mass accretion history of each halo with 


M{z) — Mo X exp [—a z] , (9) 


with Mo = M 2 oo{z — 0) and a — 2acoii- Extensions to this sim¬ 
ple function have been proposed by e. g. Tasitsiomi et al. (2004), 
McBride et al. (2009) and Correa et al. (2014), but for our purposes 
these are not necessary: the refined formulae are designed to accu¬ 
rately represent the median mass accretion histories for many halos, 
and the corrections are smaller than the scatter between individual 
halos. 

There are some differences between the conventions of Wech¬ 
sler et al. (2002) and the present work which we need to understand 
before proceeding. In the left panel of Fig. 4, the grey band and 
black dashed line show the average concentration and scatter from 
measuring the properties of ~ 120 halos of mass M ~ 10^^M© in 
our unmodified box. Wechsler et al. (2002) used a slightly different 
definition of halo mass from the one used in our study. Instead of 
defining M 200 with respect to the critical density of the universe, 
they define M^q^^ relative to the mean density. 

To show how this affects the measured concentration, the 
green points show a sample of concentration parameters estimated 
using M^o^^, using the same halos that were used to generate 
the grey shaded region. There is an overall upward offset of these 
points relative to the grey band because is correspondingly 

larger. The black solid line and the green band show the median 
relation and scatter predicted by Wechsler et al. (2002) (taken from 
their Fig. 7), corrected by a factor of 0.8 following Duffy et al. 
(2008) to account for their different erg (1, instead of 0.817 in our 
simulations). 

In summary, once the differences in conventions and cosmo¬ 
logical parameters are taken into account, we can reproduce the 
median results of Wechsler et al. (2002) in our unmodified boxes. 
The scatter from our simulation is also compatible with the much 
larger Wechsler et al. (2002) sample; outliers are likely due to the 
fact that we do not pre-select relaxed halos as they do. For consis¬ 
tency with the rest of the paper, we will use all quantities derived 
w. r. t. critical density in the following analysis. 

We are now ready to see how this relationship emerges when 
using GM halos instead of a statistical sample. For this study, we 
have selected three different halos in the reference run, which are 
all of similar mass (M ~ lO^^M©). The right panel of Fig. 4 
shows the result for 13 simulations (4 each for halos 24 & 37, and 
5 for halo 40) with different collapse times. We call each set a ‘halo 
family’, including the reference run. The right panel of Fig. 4 shows 
the results for the three families illustrated by red diamonds, green 
circles and blue squares for family 24, 37 and 40 respectively. 

The slopes of the three families appear to be consistent with, 
but scattered around, the population average (grey band). We find 
that each halo family is well-described by a linear relation 

const 1 

c =-h const 2 , (10) 

Ucoll 

which contains the offset as an additional parameter compared to 
the fit used in Wechsler et al. (2002). These fits are shown in the 
right panel of Fig. 4 as coloured solid lines; we also fit all 13 points 
together (black dashed line). In addition, the grey band shows the 
scatter expected for halos in our selected mass bin, as in the left 
panel. The fit to all 13 simulations is very similar to the median 
relation from our unconstrained box which is consistent with the 
larger sample from Wechsler et al. (2002). 

The shift of each family member along its line is dictated by 
the direction and amplitude of the density constraint in the initial 
conditions. Higher values of the density in the inner region shift 
a point towards the top left, and lower values to the bottom right 
w. r. t. the reference run. This gives us considerable insight into the 
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Figure 4. Halo concentration parameter as a function of the collapse time. Left panel: Our reference simulation gives a volume to probe the relationship using 
the traditional statistical technique. Taking 120 halos from our simulation, this results in a scatter of points in the region of the grey band. To compare with 
existing literature, we need to redefine r 2 oo (and hence c) relative to the mean (rather than critical) density, after which these halos are represented by the 
green points with error bars. The black solid line and green band show the average relation and scatter as predicted by Wechsler et al. (2002), multiplied by 
a factor of 0.8 to account for their different choice of erg (Duffy et al. 2008). Right panel: Our three constrained halos (24, 37 and 40), showing fits to each 
constrained family individually (colours) and all of them together (black dashed). The grey band is the same as in the left panel. Together the panels establish 
that (left) halos in our reference volume recover the known relationship between concentration and collapse scale factor; and (right) relationships consistent 
with this relation are also recovered individually by each GM family. The scatter of slopes between different families is expected (see text). 


kind of results that can be expected from GM compared to large 
population studies. The scatter of individual simulations within a 
GM family is very small — in other words, the concentration is 
highly predictable from a single variable. This is because, as we 
have previously emphasized, the history of each halo within a sin¬ 
gle family is as similar as possible to all the others. The normal 
scatter in the concentration-collapse relation is then seen to be due 
to factors that are not being constrained within a single family (such 
as more detailed aspects of the merger history or other variables 
such as halo spin). The GM technique allows for a detailed explo¬ 
ration of results from specific, precise changes. 

For halo 24 (red diamonds), two of the results are nearly iden¬ 
tical: the point with the highest concentration value is actually two 
points nearly on top of each other. These points have been obtained 
using the two methods for setting mass accretion history constraints 
discussed previously, emphasising that they can lead to very similar 
results. Indeed, for each halo we have performed constrained runs 
using both methods, and a mixture of the resulting datapoints are 
shown in Fig. 7. For a list of all the simulations used in this paper, 
see Table 1. 

The results of four additional runs (one each for halos 24 and 
37, two for halo 40) are excluded from this Figure. In each case, 
the estimated density profile was not well-described by an NFW 
profile due to the halo undergoing a merger or the presence of large 
substructures. This is in agreement with Zhao et al. (2003) who find 
that at least part of the scatter around the Wechsler et al. (2002) re¬ 
lation is due to poor fits to the NFW profiles and the mass accretion 
history. 


6 LIKELIHOOD OF THE MODIFIED FIELD 

As explained previously, the HR91 algorithm constructs a con¬ 
strained realization which is equivalent to (but much more efficient 
than) rejection sampling, i. e. repeatedly drawing from an ensemble 
of ACDM universes until one obtains a realization that satisfies a 
given number of constraints. However, a naive choice of constraints 
can easily result in extreme configurations which are very unlikely 
to occur within the Hubble volume of the real universe. Depending 
on context, this could even be intentional (e. g. when investigat¬ 
ing rare objects, Romano-Diaz et al. 2011a,b; Dubois et al. 2012; 
Romano-Dfaz et al. 2014); but nevertheless it is important to un¬ 
derstand how likely it is for a given constrained configuration to 
arise, relative to the reference realization. We now derive a general 
expression for evaluating this likelihood and show how it is related 
to the abundance of halos when changing the mass. 

We can compare the unconstrained and constrained fields with 
respect to the unmodified ACDM covariance matrix Co by evaluat¬ 
ing the change in defined as 

Ax" = <5n^Co“^<5„-5o+Co^5o, (11) 

where 5n is a field with n constraints. This constrained field has 
a relative abundance in the universe of compared to the 

original, unconstrained field ^o- Since this is only a relative abun¬ 
dance, applying a constraint to a halo that is rare in the reference 
simulation will in general also generate a rare object in the con¬ 
strained run. We therefore modify several halos in a similar way, in 
order to reduce the impact that picking a rare object may have on 
any of our results. 

One can calculate Ax^ directly from the density field, but we 
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Figure 5. Left panel: The relationship between Ax^ and the initial overdensity for different halo families (H24-mass, H37-mass and H40-mass; see Table 
1). Lines show the theoretical prediction from Eq. (13), whereas points give the actual change measured from the IC generator output, confirming that the 
algorithm is operating as expected. Right panel: Ax^ values (points) can be interpreted as giving the relative abundance of the halos within each genetically- 
modified family, and therefore should agree with estimates from a halo mass function (lines). The agreement is indeed good except for the H37-mass-1.5 point 
which appears to have too small a mass at z = 0 compared to expectations. This is because the halo mass function is based on an average mass build-up rate, 
whereas in this specific case the mass will only be acquired after a major merger in around 3 Gyr (see Fig. 6 and discussion in text). 


can also expand the above equation analytically by inserting Eq. (5) 
and making use of aJCoctj = 5ij. This leads to a series of cancel¬ 
lations, with the final result 

Ax^ = ^ (|di|^ - |dioO , (12) 

i=l 

where we again use dio = to express the value of the con¬ 

straints in the underlying realization. 

Crucially, the details of the original realization have dis¬ 
appeared except in the initial values of the constrained quantities, 
dio. In other words, the relative likelihood of the constrained simu¬ 
lation compared to the unconstrained case is dependent only on the 
choice of constraints. It is therefore specifically related to proper¬ 
ties of the individual halo, not to details of its surroundings. This is 
another very desirable property of the HR91 formalism and reflects 
the minimality of the changes made to the field going from to 
Sn- 

For a single constraint, Eq. (12) has a particularly transparent 
interpretation. Because of the normalization condition, the variance 
of dio for z = n = 1 in unconstrained realizations is 

(dodo) = (^Jaiaj^o) = ajCoai = 1. (13) 

Thus, for a single constraint, a change in Ax^ of 1 corresponds to 
a Icr variation in the property measured in the population-at-large. 

The left panel of Fig. 5 shows Ax^ for a single constraint as 
a function of dcons/o^ref, the ratio between a halo’s average density 
contrast after the constraint and its value in the reference run (see 
Sec. 4). The values are calculated directly from the fields (points) 
and using Eq. (12) (lines); the two methods agree to within numer¬ 
ical accuracy, which is a useful verification of the algorithm. As 


before, the results for the three families are illustrated by red di¬ 
amonds, green circles and blue squares for haloes 24, 37 and 40 
respectively. The minimum at dcons = 0, as well as the symmetry, 
is expected for a zero-mean Gaussian random field. 


6.1 Connecting initial conditions and non-linear structure 


In this section, we work towards establishing a quantitative connec¬ 
tion between the degree of change in the initial conditions and in 
the final, non-linear structure. Using a single constraint, we inves¬ 
tigate how a change in density contrast in the initial conditions is 
related to the resulting halo mass at late times. Qualitatively, an in¬ 
crease in overdensity should lead to a more massive object at z = 0 
as explained in Sec. 4. 

Quantitatively, the probability of finding a halo with mass M 
at z = 0 is given by the halo mass function, n(M)dM, which de¬ 
pends on the cosmological power spectrum and growth function. 
Given two halos of mass Mo and Mi, their relative abundance 
is given by the ratio of the halo mass function at those masses, 
n(Mi)/n(Mo). Assuming that the statistical properties of con¬ 
strained and unconstrained simulations can be related by the change 
in mass of the target halo alone, this ratio also gives the relative 
probability of the structure in the two simulations. 

Additionally, we can calculate the relative probability of the 
two initial conditions using Eq. (11); specifically 


P{di) 

p{do) 


exp [-Ax^/2] , 


(14) 


where p{d) is short-hand for the probability of a constrained field 
using one ‘density constraint’ with value d. Since we only consider 
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Figure 6. Slices from the original simulation (upper panel) and H37-mass- 
1.5 simulation (lower panel) illustrate how the target halo is, at 2 : = 0, seen 
at a time where it is about to undergo a major merger in the latter case. 
For that reason its mass undershoots the expectation from the halo mass 
function (Fig. 5) which averages over all the possible discrete realizations 
of the accretion history. Black circles show the size of the virial radius. 


probability ratios and Ax^, any terms which only change the nor¬ 
malization of p{d) have dropped out. 

Now we have two methods of calculating the relative proba¬ 
bilities, and we can check whether they agree. We rewrite p{d) in 
terms of halo mass using the conservation of probability 

p{M)^p(d)d'{M), (15) 

where M is the halo mass at z = 0, and d'{M) is the derivative 
of the constraint d with respect to the mass, evaluated at M. By 
combining Eqs. (14) and (15), we can find a relationship between 
thop{Mi)/p{Mo) and Ax^; namely 

= exp [-AxV2] d'(Mi) [d'(Mo)]“'. (16) 

If our assumptions are correct, this expression for p{Mi)/p{Mo) 
should be equal to the mass function ratio n{Mi)/ n(Mo). We have 
used HMFcalc (Murray et al. 2013) to generate a halo mass func¬ 
tion at z = 0 for our cosmological model, and confirmed that it 
provides a good fit to our simulations.^ 

Evaluating Eq. (16) also requires an estimate of the Jacobian 
factors on the right hand side. This can be obtained either from a 
physical model underlying the mass function or by using an empiri¬ 
cally calibrated M{d) relationship from the simulations. We chose 
the latter approach by fitting a power-law relation between d and 
M, which allows us to obtain values for d'{M) at different halo 
masses separately for each halo family (24, 37 & 40; introduced in 
the previous section). This leaves us with a ‘semi-analyticaT pre¬ 
diction: theoretical halo mass function plus fit to the Jacobian. 


The right panel of Fig. 5 shows the results of the calcula¬ 
tion (lines), as well as points evaluated directly from the simula¬ 
tions. Overall, most points show a good agreement: there is con¬ 
sistency between the population statistics and the abundance calcu¬ 
lated from the GM Ax^ values. The broad agreement justifies our 
set of assumptions for calculating abundances in this specific case 
of a single density constraint. However the individual halos do scat¬ 
ter around the relation and there is one point that clearly does not 
fit the expectations. This arises from our H37-mass-1.5 run where 
we increase the initial overdensity of the proto-halo 37 region by a 
factor of 1.5. 

The mismatch can be understood by considering the discrete 
nature of merger histories. Specifically, Fig. 6 shows the projected 
density at the last output (z = 0) in a region around halo 37 in the 
original run (upper panel) and the xl.5 run (lower panel). In the 
latter case, a major merger (mass ratio ~ 2) will occur in around 
3 Gyr. After this merger, the anomalous point will shift signifi¬ 
cantly rightwards in Fig. 5 to the correct mass ratio 3.1) ac¬ 
cording to the Ax^(M) derived from the halo mass function.^ We 
can frame this in another, more general way: the halo mass func¬ 
tion is a statistical construction that corresponds to averaging over 
all possible histories, but the individual points in Fig. 5 represent 
modifications to specific halos which have a discretized accretion 
history. Therefore, they scatter away from the line, especially when 
seen at special times (such as shortly before a major merger). 

The main conclusion from Fig. 5 is therefore that the changes 
in the x^ give us a good quantitative handle on the relative abun¬ 
dance of halos of different types, at least in this case where we have 
only changed the mass. However, the complex non-linear connec¬ 
tion between initial and final states means that Ax^ will always 
need to be interpreted with care. 


7 DISCUSSION & CONCLUSIONS 

In this paper we have demonstrated an extension of the HR91 tech¬ 
nique to modify the initial conditions of a numerical simulation. 
For this modification, we selected regions of arbitrary shape, de¬ 
fined solely by the particles that form a halo in our reference simu¬ 
lation. This is a different approach than that used in previous works, 
which relied on imposing constraints of a given analytic profile. By 
applying our constraints only to the halo particles, we showed that 
we can ‘genetically modify’ a single object, changing its properties 
in a smooth and continuous way. 

Using constraints on the density averaged over all halo parti¬ 
cles controls the total mass, whereas adding additional constraints 
allows us to change the halo’s collapse time. This serves as a 
demonstration of the technique and is the basis of the creation of 
further constraint types to study the impact of other halo character¬ 
istics on a halo’s evolution. 

Using the collapse time constraint, we investigated the den¬ 
sity profiles of the resulting halos at z = 0 and found that the 
distribution of their concentration parameter is consistent with the 
results of statistical analyses such as Wechsler et al. (2002). How¬ 
ever, we also find that different halos occupy different regions in the 
parameter-space, and have different trajectories when their collapse 
time is changed. We plan to study this behaviour in future work, 
in order to determine which other halo parameters have changed. 


^ This is preferable to obtaining a (noisy) estimate of the mass function 
directly from our limited volume; the fitting functions included in the HM¬ 
Fcalc tool were validated using detailed studies of large simulation suites 
(e. g. Tinker et al. 2008). 


^ Note that the nearest massive halo in the original run corresponds to the 
same particles, but is considerably further away from halo 37 and is not on 
a trajectory that will lead to a merger within a Hubble time. 
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This should be complementary to the principal component analy¬ 
sis carried out by Skibba & Maccio (2011); Jeeson-Daniel et al. 
(2011); Wong & Taylor (2012), which revealed somewhat incon¬ 
clusive correlations between additional internal halo parameters. 
With our constrained simulations, we will not only be able to find 
correlations but to explicitly test their significance. Since we can 
directly compare the constrained halo to its reference in the uncon¬ 
strained run, we can establish exactly which changes in the halo 
parameters have a physical impact. 

We have provided a way of quantifying the likelihood of mod¬ 
ified initial conditions via a difference between the constrained 
and unconstrained fields. In general, this statistic can be used to 
assess how compatible the modified object is with the underlying 
cosmology. Similar expressions were obtained by van de Weygaert 
& Bertschinger (1996); however our orthonormalisation procedure 
allows for the derivation of the considerably simpler Eq. (12). The 
statistic can be used to quantify the rarity of genetically modified 
objects relative to the unconstrained realization. As an example we 
showed that modifying the mass produces abundance constraints 
that are quantitatively consistent with the traditionally-measured 
halo mass function at z = 0. Individual halos have discrete ac¬ 
cretion histories and scatter around the mean relation predicted by 
the mass function; the strongest outlier in our study is about to un¬ 
dergo a merger at z = 0, which significantly lowers its current 
mass. In principle, the Ax^ measure could also be used to specifi¬ 
cally create objects that are ‘rare’ in a ACDM universe (similar to 
e. g. Romano-Diaz et al. 201 la,b; Dubois et al. 2012; Romano-Diaz 
et al. 2014); we leave such a study to future work. 

In this paper we have used a uniform resolution over a box 
size of 50 h~^ Mpc; having a sufficiently large box is important 
to ensure that halos are embedded in the correct large-scale en¬ 
vironment. Our code also produces ICs for ‘zoom’ simulations 
with varying resolution (see also Prunet et al. 2008; Romano- 
Diaz et al. 2014); the only major difference when generating these 
is the extra computational complexity introduced in the transfor¬ 
mation between real space and Fourier space on an irregularly- 
spaced grid, which has been tackled elsewhere in the literature (e.g. 
Bertschinger 2001; Hahn & Abel 2011). The real power of the ap¬ 
proach to generate insight into a population from a handful of runs 
will become more apparent as we begin to use these zoom ICs in 
tandem with high-resolution baryonic physics. 

While our focus here has been on basic properties such as the 
formation time and mass of a system, many other interesting as¬ 
pects of evolution can be changed by constraining different proper¬ 
ties. One example with which we are experimenting is the specific 
angular momentum, which can be controlled because tidal torque 
theory describes the connection between ICs and final spin (e.g. 
White 1984; Catelan & Theuns 1996; Porciani et al. 2002a,b); the 
internal properties of the galaxy forming inside the dark matter halo 
will naturally depend on the spin parameter of the halo. We are able 
to generate initial conditions that modify the spin parameter of the 
halo, but leave the mass and merger history untouched. Studies of 
halo spin constraints for dark matter and hydrodynamic simulations 
will be presented in future work. 
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APPENDIX A: TECHNIQUE 

Here, we will give a detailed description of the derivation of the 
HR91 operator and our Eq. (5). 

The values of the three-dimensional overdensity field 6 = 
S{x) in the initial conditions of our cosmological simulations are 
distributed according to a multivariate Gaussian 

po{S) oc exp - Mo)^Co ^(<5 - j , (Al) 

with mean {S) = fj,Q and Co = {{S — {S — /ito)) (the covari¬ 
ance, or the power spectrum in Eourier space). Cosmological initial 
conditions have zero mean, /k-q = 0, but we will consider the fully 
general case. 

We will build the general procedure by induction. Suppose we 
have SL Pi -1 (5) that describes the probability distribution function 
for i—1 constraints; we now want to add the Ah constraint, ensuring 
that Ck-ld = di for some constraint vector cXi and constant di. To 
gain samples from the constrained distribution one could sample 
from the original distribution and reject all those trials which lie 
I t P 

too far away from — dJ = 0. Mathematically this can be 
expressed by multiplying the original probability distribution by a 
penalty function, e. g. 

Pi{S) (X lim pi_i(^) exp ( laj^ - dJ ), (A2) 

/0^oo y z I \ J 

where the constant of proportionality renormalizes the probability 
distribution function and is dependent on jd. In the limit /S ^ oo, 
the penalty function becomes a Dirac-Delta distribution and the 
constraint is satisfied exactly. 

Under the assumption that pi-i is Gaussian, the new prob¬ 
ability function is the product of two Gaussians, and so remains 
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Gaussian itself; consequently after imposing i constraints we must 
be able to write 

Pi{5) oc exp (A3) 


for some mean ^jl- and covariance Ci which we will now derive. By 
multiplying out Eq. (A2) we obtain 


Pi{S) (X lim exp 

j3^oo 

~ fA + 


(c-_\ + ^3cx^cx]) {S - 
^ , (A4) 


where we have already thrown away several terms which are zero- 
order in 5 since they just change the normalization. By compar¬ 
ing terms in Eqs. (A3) and (A4) we can first read off = 
C~\+^OLiOL \. We will also need a normalization for the a^, which 
conveniently can be chosen^ as 

cx\Ci-icxi — 1. (A5) 


Next, we apply the Sherman-Morrison formula (Sherman & Mor¬ 
rison 1950), 

/f— — 1 I n t\—1 f” Q —1 

( Q _1 + ^oliol\) = Q -1 - / 3 --- 

1 + l3cx\Ci-icXi 




(A6) 


where we have used /3 ^ 1 and the normalization condition (A5) 
in the second step. 

The terms in the second line of Eq. (A4) have to cancel exactly. 
Plugging Eq. (A6) into this expression leads to 

Mi = lim Mi-1 - (1 - /3“^)Ci-iaialMi-i + dia.iCi-iot\ 

p^oo 

= Mi-1 + Ci-ifti {di - alMi-i) > (A7) 

and finally taking the limit in Eq. (A6) yields 

Ci = Ci—1 Ci—iQ^iQiJCi—1. (A8) 

This result allows us to apply as many constraints as desired analyt¬ 
ically - by looping over the constraints and updating the covariance 
matrix and mean at each step, then drawing a constrained realiza¬ 
tion - but this would be computationally expensive. Instead, the 
constrained realization can be constructed from the unconstrained 
field using a projection operator, which we will now derive. 

Eor notational simplicity, in addition to normalizing the con¬ 
straints, it is also helpful to make them orthogonal (e. g. through 
a Gram-Schmidt procedure) in the sense that aJCoctj = 0 for 
i ^ j. One can then verify by substitution (see Sec. Al) that the 
constrained field has mean 

n 

Mn = Mo + Coa, (di - a^Mo) (A9) 

i=l 

and covariance 

n 

C„ = Co - Coaia+Co (AlO) 

i=l 

for orthonormalized {a^}. 

Efficiently drawing from the distribution implied by the above 


Unless OLi is a null direction of Cj-i, but then there would be zero prob¬ 
ability of our constraint in the original distribution. 


mean and covariance is made possible by any operator Pn that takes 
a realization from the unconstrained field and forms a new real¬ 
ization under n constraints via the ansatz 

Sn = Pn (^0 - Mo) + Mn’ (^H) 

where to gain the correct covariance Cn = ((^n — Mn) — Mn) ^) 

one must demand 

P^CoPt = Cn. (A12) 

There are an infinity of operators Pn with this property: Given 
any specific Pn one can form P^ = UPn where U^CqU = 1, 
and the new P^ satisfies the required identity (A 12). To obtain the 
unique HR91 operator, we additionally require Pn to make min¬ 
imal changes to the field. This implies Pn^n = ^n — in other 
words, that no changes are made if the field already satisfies the 
constraints. Using Eq. (All), it immediately follows that P^ = Pn 
and PnMn — Mn- of these conditions implies that all 

eigenvalues of Pn are either 1 or 0. 

One can verify by substitution that all these requirements are 
satisfied by 

n 

Pn = 1 — CoCXiOil (for orthonormalized {ai}), (A13) 

i=l 

Note that the HR91 form given in their Eqs. (2) - (4) builds the 
orthonormalization procedure into the projection operator (appear¬ 
ing as in their notation). However, as stated above we found it 
notationally simpler to pre-condition the constraints into orthonor¬ 
mal form using the Gram-Schmidt procedure. Both formulations 
are mathematically equivalent (see Appendix B). 

Inserting Eqs. (A 13) and (A9) into (All) then leads to the final 
expression 


n 

= ^0 + ^ CoCKi {di - dio ), (A14) 

i=l 

where dio = 

In practice, most of the necessary calculations are performed 
in Eourier-space, because there Co is the ACDM power spectrum 
which is diagonal. Any constraint vector cti and density field S can 
be easily converted using numerical East Eourier transformations. 

Note that the algorithm in its current form only takes into ac¬ 
count the contribution from the power spectrum. If one wanted to 
generate constrained initial conditions based on an observational 
dataset, the associated uncertainties would introduce extra contri¬ 
butions in the new covariance matrix (e. g. Zaroubi et al. 1995; van 
de Weygaert & Bertschinger 1996), which is not included in the 
current implementation. 


Al Comments on normalisation 

Throughout this paper, we use the same notation for the normalised 
and unnormalised constraints (expressed by = di). In prac¬ 
tice, these quantities are affected by the normalisation condition in 
the following way: if aJCocti = hn before normalisation, then we 
immediately find cxi cxils/^i to satisfy aJCoCKi = 1. Accord¬ 
ingly, the constant di transforms as di di /as well, such 
that cxld = di is still obeyed after normalisation. 

The consistency of the Gram-Schmidt condition aJCoCKj = Sij 
and our normalisation aJCi-ia^ = 1 can be shown as follows: 
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Multiplying (A 10) by an on both sides yields 
^n^n — lOln — aJj^Coan 

n —1 

- ^ aJ^Coa^aJCoan, (A15) 

i=l 

where we have shifted the index n by 1 eompared to (A 10) for ease 
of notation. Inserting ctj Co cxj = Sij into the right hand side leads 
to 

n —1 

^n^n — lOLn — 1 ^ ^ ^ni^in — I 7 (A16) 

i=l 

which proves the equality. 


APPENDIX B: EQUIVALENCE TO HR91 

In this section we translate our notation into the one used by HR91 
to show that our final expression (A14) is equivalent to their Eq. (4), 
which states 

n 

Sn = So+J2Ul'/idi-dio), (Bl) 


where we have already translated their notation for the constrained 
and unconstrained field, and the values of the constraints (f ^ S, 
Ci di). The two additional functions are 


= (<5oCt> 

(B2) 

and 


in = {CiC]), 

(B3) 

where © = aj^o (not to be confused with our covariance ma¬ 
trix Ci), and we have added daggers to the HR91 notation to al¬ 
low for complex-valued quantities. Inserting these expressions into 
(Bl) leads to 

—1 

= ^0 + ^ Coai l^aJCoa^j {di - 

dio) (B4) 

n 

= ^0 + ^ ^ CoOLi{di — dio\ 

(B5) 


i=l 


where in the first step we used that Co = (^o^o) because /itg = 0 
in HR91, and our orthonormalisation condition in the second step. 
This shows that the Gram-Schmidt approach is equivalent to per¬ 
forming the matrix inversion ^. 
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